Dynamics and morphology of dendritic flux avalanches in superconducting fllms 
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We develop a fast numerical procedure for analysis of nonlinear and nonlocal electrodynamics of 
type-II superconducting films in transverse magnetic fields coupled with heat diffusion. Using this 
procedure we explore stability of such films with respect to dendritic flux avalanches. The calculated 
flux patterns are very close to experimental magneto-optical images of MgB2 and other supercon- 
ductors, where the avalanche sizes and their morphology change dramatically with temperature. 
Moreover, we find the values of a threshold magnetic field which agrees with both experiments and 
linear stability analysis. The simulations predict the temperature rise during an avalanche, where 
for a short time T ^ 1.5Tc, and a precursor stage with large thermal fluctuations. 

PACS numbers: 74.25.Ha, 68.60.Dv, 74.78.-w 



I. INTRODUCTION 

The gradual penetration of magnetic flux in type- 
II superconductors subjected to an increasing applied 
field or electrical current can be interrupted by dramatic 
avalanches in the vortex matter P The mechanism respon- 
sible for the avalanches is that an initial fluctuation re- 
duces locally the pinning of some vortices, which start to 
move, thus creating dissipation followed by depinning of 
even more vortices. A positive feedback loop is formed 
where a small perturbation can escalate into a macro- 
scopic thermomagnetic breakdown.^ 

In thin film superconductors, the dynamics and mor- 
phology of these avalanches is tantalizing, when at 
very high speeds they develop into complex dendritic 
structures, which once formed remain robust against 
changes in external conditions. When repeating iden- 
tical experiments one finds that the patterns are never 
the same although qualitative features of the mor- 
phology, such as the degree of branching and overall 
size of the structure, show systematic dependences on, 
e.g., temperature. Using magneto-optical imaging flux 
avalanches with these characteristics have been observed 
in films of Nb,^ YBasCuaOr-^J MgBs,"^ NbsSn^ 
YNi2B2C,'^ and NbN.^^ Investigations of onset conditions 
for the avalanche activity have identified material depen- 
dent Uireshold values in temperature,^ applied magnetic 
field, and transport current f^l as well as in sample 
size.^^ Analytical modeling of the nucleation stage has 
explained many of these thresholds using linear stability 
analysis .1^2^ 

Far from being understood is the development of the 
instability from its nucleation stage to the fully devel- 
oped dendritic pattern. Aranson et al.^^ explored the 
dynamics of the flux avalanches as a numerical solution 
of Maxwell's equations with temperature dependent crit- 
ical current density. The dynamical process was gov- 
erned by the interplay between an extremely nonlinear 
current-voltage relation, heat diffusion, and the nonlocal 
electrodynamics characteristic for thin superconducting 
films. To treat the nonlocal electrodynamics the authors 



used periodic continuation of the sample taken as an infi- 
nite strip. This scheme should be a good approximation 
inside the sample, although not necessarily close to the 
edges. In fact, in thin films the magnetic field near the 
edges is significantly enhancecff^ due to the flux expul- 
sion. Moreover, all experiments show that the instability 
is always nucleated at an edge. Therefore, a careful ac- 
count of the electrodynamics close to the edges, including 
the regions outside the film, is expected to be crucially 
imp or t ant 

In this work we study the formation and characteris- 
tics of dendritic flux avalanches using a numerical scheme 
that takes into account the nonlocal electrodynamics 
both inside and outside a finite-sized superconducting 
film. It is shown that our simulations largely reproduces 
experimental results obtained by magneto-optical imag- 
ing of dendritic avalanches in films of MgB2, and further- 
more gives detailed insight into not yet observed quanti- 
ties such as local temperature rise and electrical field. 

The paper is organized as follows. Section [ll| presents 
the model and the equations describing the process. 
The numerical scheme including the implementation of 
boundary conditions and thermomagnetic feedback is de- 



scribed in Sec. |III| The results for the time-dependent 
distributions of magnetic flux and temperature are pre- 
sented and discussed in Sec. |IV[ while Sec. [V| gives the 
conclusions. 



II. MODEL 

Consider a rectangular superconducting film zero-field 
cooled below the critical temperature, Tc, followed by 
a gradual increase in a perpendicular applied magnetic 
field. The film is deposited on a substrate, which in the 
process will be regarded as a sink for the dissipated heat. 
Shown in Fig. [l] is a sketch of the overall configuration, 
including the relevant fields and currents. 

The macroscopic behavior of type-II superconductor 
films in a transverse applied magnetic field. Ha, is well 
described by quasi-static classical electro dynamics .^^^^ 
Here the sharp depinning of vortices under flowing cur- 
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FIG. 1. (Color online) Schematic of the sample configuration. 



rent is represented by a highly nonlinear current-voltage 
relation 



E = p(J)J/d, 



P{J) 




J <Jc, T <Te, 
J>Jc, T<Te, 
T>Tr,. 



(1) 



Here E is the electric field, J is the sheet current (J = 
|J|), Jc the critical sheet current, n is the creep exponent, 
po is a resistivity constant, pn is the normal resitivity, 
and T is temperature. It is assumed that the sample 
thickness, d, is so small that variations in all relevant 
quantities across the thickness can be ignored. For T < 
Tc the temperature dependence of the critical current and 
flux creep exponent^^ are taken as 

Jc = Jco{l-T/T,) and n - 1 = rioTjT , (2) 

where Jco and no are constants. 

The distribution of temperature is described by the 
heat diffusion equation 



dcf = dV ■ [kS/T) - /i(T - To) + J • E , 



(3) 



where n is the thermal conductivity of the superconduc- 
tor, c is its specific heat. To is the substrate temperature, 
taken to be constant, and h is the coefficient of heat 
transfer between the film and the substrate. The k.^ c 
and h are all assumed to be proportional to T^, whereas 
a relatively weak t emperature dependences of po and pn 
are neglect ed.'^^'^ 

Following Ref. 17 we define the local magnetization, 
g g{r), as 



X z = V X (gz) 



(4) 



where r = (x^y) is a 2D vector in the film plane, and z is 
the unit vector in the perpendicular direction. Outside 
the sample there are no currents, and we set ^ = by 
definition. The Biot-Savart law can then be written as 



Mo 



Ha = Qg= / d^r'Q{r-r',z)g{r'). 



(5) 



where the integral is calculated over the whole plane. The 
kernel Q{t) should be calculated as a limit at z ^ of 
the expression 



Q{r,z) 



1 2^2 



47r (^2 + ^2)5/2 



(6) 



Here reqularization is needed to avoid formal divergence 
of the r.h.s. of Eq. ([5| at 2: = 0, r = r^ The Fourier 
transform of lim^^o z) is equal to Therefore, 
from the convolution theorem it follows that the inverse 
operator acting on some function (p{r) can be ex- 
pressed as 



(7) 



Here F[(^{y)] and [(/:?( k)] are Fourier and inverse 
Fourier transform, respectively, and k = |k|. 

Inverting Eq. ([5| one arrives at the equation for the 
time evolution of the local magnetization, 

^(r,t) = 2^-i(fc-i^[/io"'^.(r,t)-i^a(t)]) . (8) 

Equations ([3|, Q and ([8| therefore determine the dy- 
namics of ^(r,t), T(r, t), etc. To solve these equations 
numerically we proceed from the continuous to a discrete 
formulation. 



III. NUMERICAL APPROACH 

To allow use of the fast Fourier transform (FFT) we 
consider a rectangular area of size 2Lx x 2T^ contain- 
ing the sample plus a substantial part of its surrounding 
area. A key point is to select proper values for and Ly 
relative to the sample size, 2a x 2h. By including too little 
area outside the sample one clips away the slowly decay- 
ing tail of the stray fields, leading to decreased accuracy 
at large scales, and major deviations from the correct 
physical behavior.^ On the other hand, including too 
much of the outside area keeping the same number of 
the grid points tends to decrease the accuracy at small 
scales, where actually the most interesting features of the 
dendritic avalanches appear. This blurring can be com- 
pensated by using a finer spatial grid, at the cost of a 
rapidly increasing computation time. 

A careful test of our numerical scheme was done by 
comparing the calculations with the exact solution for 
the Bean critical state in an infinitely long strip. It is 
found that already with Lx/a > 1.3 the calculated re- 
sults are correct within a few percent, and are essentially 
indistinguishable from the exact solution in graphic com- 
parisons. 

In the FFT-based calculations the rectangle 2Lx x 2Ly 
is discretized as a A^^^ x Ny equidistant grid, and used 
as unit cell in an infinite super lattice. The Fourier 
wave vectors kx^y are then discrete, kx^y = '^(Ix.yl ^x.y, 
where qx^y are integers. The Brillouin zone is chosen as 
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\(lx,y\ ^ which ensures g{r^t), T{T^t), etc. to be 

real valued. 

The calculation of the temporal evolution is based on 
a discrete integration forward in tim^^ of the local mag- 
netization 



(9) 



starting from ^(r, 0) = 0. Once giv^t) is known at time 
t, we proceed one time step by determining g{r,t). The 
g{r^t) can be calculated from Eq. ([8|, provided is 
known everywhere within the unit cell. For this, we have 
to find self-consistent solutions for g and given the 
function g. 

For the area inside the superconductor the material 
law, Eq. ([T]), applies and together with the Faraday law, 
= -(V X E)^, it follows that 



B, = \/-{p\/g)/d. 



(10) 



The gradient Vgir^t) is readily calculated, and since the 
result allows finding J(r, t), from Eq. Q, also p{r^t) is 
determined from Eq. ([T]). The difficult point is that g 
depends on the distribution of B^ in the whole unit cell. 
The task is to find the Bz outside the sample which leads 
to ^ = outside. This cannot be calculated directly since 
there is a nonlocal relation between Bz and g. Instead 
we use an iterative procedure. 

Let us label the iterations by a superscript (i). At the 
first step, i = 1, we calculate Bz inside the superconduc- 
tor from Eq. (10). Then an initial guess is made for the 

time derivative, Bi^\ outside the sample. From Eq. (|8| 
we now compute the time derivative g^^\ In general, 
this ^f*^^^ does not vanish outside the superconductor. To 
correct for this, a new and improved Bz is chosen as 

Bi'^^^ = B^^ - iJLoQOg^^ + C^^, (11) 

where the projection operator O vanishes inside the su- 
perconductor and equals to 1 outside it. The constant 
C^*^ is determined by the fiux conservation. 



(12) 



The procedure is stopped after s iterations when the val- 
ues of g outside the superconductor becomes sufficiently 
small. The final distribution, ^^*^(r), is taken as the 
"true" g{r,t)^ and substituted into Eq. ([9| in order to 
advance in time. 

A good choice for the initial state of the iteration at 
time t is Bi^\t) = Bi^\t — At)^ i.e., each iteration starts 
from the final distributions achieved during the previous 
iteration. Normally, s = 5 iterations is sufficient to give 
good results. 



IV. RESULTS AND DISCUSSION 

Numerical simulations were performed for samples 
shaped as a square of side 2a and with an outside area 




FIG. 2. (Color online) Calculated distribution of Bz at 
an applied field Ha = O.lSJco, and substrate temperature 
To = Tc/4:. The image brightness represents the magnitude of 
Bz . The sample contour appears as a bright rim of enhanced 
field, and the black central area is the flux-free Meissner state 
region. 



corresponding to = Ly = 1.3a. The total area is dis- 
cretized on a 512x512 equidistant grid. Quenched disor- 
der is included in the model by a 10% reduction of Jco 
at randomly selected 5% of the grid points. The simu- 
lated ffux penetration process starts at zero applied field 
with no fiux trapped in the sample, which has a uniform 
temperature Tq. 

Calculations were performed at Tq = Tc/4 using 
material parameters corresponding to a typical MgB2 
film,C^Pn=7 fincm, n = 0.17 kW/Kmx (T/TJ^ and 
c = 35 kJ/Km^ x (T/Tc)^, where pn is the normal resis- 
tivity at Tc = 39 K, Jco = 50 kA/m, po = p^, d = 0.5 //m, 
a 2.2 mm, and h 220 kW/Km^ x {T/Tcf. We choose 
no = 19 and limit the creep exponent to n(T) < n^ax = 
59. The field was ramped from Ha = at a constant 
rate. Ha = lO'^JcoPn/adpo. 

Figure [2] shows the ^^-distribution at poHa = 
O.lSpoJco = 11 mT, where three large dendritic struc- 
tures have already been formed. The numerical labels 
indicate the order in which they appeared during the field 
ramp. The first event took place at the threshold applied 
field, poHth = 0.145/io^cO = 9.1 mT, which is in excellent 
agreement with measurements on MgB2 films just below 
10 Tc/4:. At lower fields, the fiux penetration was 
gradual and smooth, just as seen on the left edge of the 
sample, where the characteristic "pillow effect" for films 
in the critical state is very well reproduced. 



FIG. 3. (Color online) Map of the sheet current, J, corre- 
sponding the image in Fig. |2] The brightness represents J, 
where black means J = 0. 



The dendritic avalanches all nucleate at the edges, and 
one by one they quickly develop into a branching struc- 
ture that extends far beyond the critical-state front and 
deep into the Meissner state area. The trees are seen to 
have a morphology that strongly resembles the flux struc- 
tures observed experimentally in many superconducting 
films .'^''^ The simulations also reproduce the experimental 
finding that once a flux tree is formed, the entire dendritic 
structure remains unchanged as Ha continues to increase. 
The supplementary material^^ includes a VIDEO clip of 
the dynamical process, and shows striking resemblance 
with magneto-optical observations of the phenomenon. 

Figure |3] shows the sheet current magnitude, J, corre- 
sponding to the flux distribution in Fig. [2] From this map 
it is clear that the dendrites completely interrupt the cur- 
rent flow in the critical state, and redirect it around the 
perimeter of the branching structure. This vast perturba- 
tion of the current has been demonstrated experimentally 
earlier using inversion of magneto-optical images. Note 
that the critical state region contains dark pixels which 
are the randomly distributed sites of reduced Jco- 

To investigate reproducibility in the pattern formation, 
microscopic fluctuations were introduced by randomly al- 
ternating between right- and left-derivatives in the dis- 
crete differentiation. Due to the nonlinear form of Eq. ([T]) 
this procedure gives large local variations in the electrical 
fleld. Figure |4] shows an overlay of two simulation runs 
with different realizations of the microscopic fluctuations 
while keeping the same quenched disorder in Jco- The 
two resulting images were colored so that adding them 



FIG. 4. (Color online) Color-coded overlay of two separate 
runs with same quenched disorder but with different micro- 
scopic fluctuations. The pixels in gray-scale represent overlap- 
ping results. The parameters are the same as in the caption 
of Fig. [I 



gives shades of gray where both coincide in pixel val- 
ues. Clearly, the two runs gave different results as far as 
the dendritic pattern is concerned. Both produced three 
branching structures, where two are rooted at the same 
place and the third is at a different location. Even for 
those with overlap, there are parts of the structure that 
differ considerably, especially in the flner branches. In 
contrast, both the critical state and the Meissner state 
regions are essentially identical in the two runs. Note 
the color at the edge of the right hand side near the root 
of the green dendrite, which reflects that the growth of 
the flux structure drains the external field near the root. 
Moreover, the root of all the trees are not far from the 
middle of the sides. Both features are in full accordance 
with experiments. 

Each dendritic avalanche is accompanied by a large lo- 
cal increase in temperature. Shown in Fig. [5^ is a plot 
of the maximum temperature in the film during a field 
ramp with substrate kept at Tq = Tc/4. The spikes in 
the temperature rise as high as 1.5Tc. The maximum 
temperature is found in the root region of the avalanche. 
The heating above Tc is an interesting prediction; to our 
knowledge, the temperature of propagating avalanches 
has not been observed experimentally. At the same time, 
the result is consistent with the measured heating of uni- 
form fiux jumps in Nb foilJ^ and the magnetic field- 
induced damage in a YBa2Cu307_a; film during dendritic 
growth.^ 
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FIG. 5. Maximum temperature in the superconductor during 
an ascending field ramp at To = Tc/4. The panels (a)-(c) are 
successive magnifications of the first avalanche event. 

The first avalanche in Fig. [5^ appears at i^th = 
0.145Jco. Since the chosen disorder is rather weak and 
the ramp rate is high, the heat diffusion to the substrate 
is expectedly a more important stabilizing factor than 
lateral heat diffusion, the theoretically predicted thresh- 
old field i&i2 

i7th = ^tanh-^f . ) . (13) 

At T = Tc/4 and with n = 59 this gives Hth = O.lSJco, 
in excellent agreement with the present simulation. Here, 
= 0.6Jco is the effective critical current, which is lower 
than Jc due to fiux creep. At the same time, the adiabatic 
threshold field^^ is much smaller than i^th, which means 
that the heat diffusion and heat transfer to the substrate 
prevent avalanches. However, during short time intervals 
cooling is not always effective, and the temperature ex- 
periences large fluctuations. The fluctuations are partic- 
ularly large as Ha approaches triggering of an avalanche, 
see Fig. ^p. In these intervals both heat absorption and 
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FIG. 6. (Color online) Magnetic moment in units of mo = 
a^Jco as function of increasing field obtained by simulations 
at three different temperatures, To. Each jump in the curves 
represents a flux avalanche. 



lateral heat diffusion play important roles in stabilizing 
the superconductor. A close-up view of the maximum 
temperature during the first avalanche at Tq = Tc/4 
is shown in Fig. [Sj^. First, the temperature rapidly in- 
creases, and then decays much slower. The duration of 
the avalanche is 0.18 /is. Since the length is 2.5 mm, the 
average propagation velocity is of order 14 km/s. This 
numerical value is reasonable compared to previous mea- 
surements, where the flux dendrites were triggered by 
a laser pulse in YBaCuO films.I^ The maximum electric 
field in the superconductor during the avalanche is also 
high, found from the simulations to be approximately 
5 kV/m. 

The abrupt redirection of the current implies that the 
magnetic moment of the sample makes a jump and be- 
comes smaller. Figure [6] shows the moment as function 
of the increasing applied field. Each vertical step corre- 
sponds to a fiux avalanche. The lower curve, obtained for 
To = Tc/4, shows jumps with typical size of O.lmo with 
a slight dispersion, which is due to variations both in 
shape and location of the avalanches. More pronounced 
is the variation in jump size with temperature. As Tq 
gets lower the jump size becomes smaller, and the events 
more frequent. In the graphs for Tq/Tc = 0.20 and 0.17, 
the jump size reduces to 0.03mo and O.Olmo, and jumps 
appear on average with field intervals of AHa/ Jco = 0.01 
and 0.002, respectively. In real samples a similar temper- 
ature variation of jumps in the m-H curves was observed 
by magnetometry.'2'^^ESI 

It has been reportecP that the morphology of flux 
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avalanches is strongly temperature dependent. This is 
illustrated in the bottom panel of Fig. [7| showing three 
magneto-optical images of a 0.4 /im thick MgB2 square 
film at To =4 K, 6.3 K and 7.9 K. The images show a 
crossover from many long fingers at 4 K to medium sized 
dendrites at 6.3 K, to a single highly branched structure 
at 7.9 K. The simulation results shown in the top panels 
reproduce this result and show exactly the same trend 
as the experiments. At the lowest temperature, 0.17Tc, 
there are many finger-like avalanches. At the middle tem- 
perature 0.2Tc there are fewer avalanches, with typically 
three to four branches each. At the highest temperature 
0.25Tc there is just one big avalanche, with seven main 
branches. 



V. CONCLUSION 

In conclusion, we have developed and demonstrated 
the use of a fast numerical scheme for simulation of 



nonlinear and nonlocal transverse magnetic dynamics 
of type-II superconducting films under realistic bound- 
ary conditions. Our simulations of thermomagnetic fiux 
avalanches qualitatively and quantitatively reproduces 
numerous experimentally observed features: the fast fiux 
dynamics, morphology of the fiux patterns, enhanced 
branching at higher temperatures, irreproducibility of 
the exact fiux patterns, preferred locations for nucleation, 
and the existence of a threshold field. The scheme allows 
determination of key characteristics of the process such 
as maximal values of the temperature and electric field 
as well as typical propagation velocity. 
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